********************************************************************************
* CRISIS WINDOW ROBUSTNESS TESTS WITH SMALL-SAMPLE CORRECTIONS
* 
* Purpose: Test sensitivity of phase-specific FAC influence estimates to 
*          alternative crisis window definitions (18, 24, 30 months)
*          with appropriate small-sample inference (HC3, bootstrap)
*
* Output: Tables B4-B6 for Appendix B
*
* Author: O'Flynn & Stockhammer
* Date: December 2024
********************************************************************************

clear all
set more off
cap log close
log using "window_robustness.log", replace

* Install required packages
cap which outreg2
if _rc ssc install outreg2, replace

********************************************************************************
* SETUP
********************************************************************************

cd "/Users/ciaran/PycharmProjects/fac"

cap mkdir "results"
cap mkdir "results/robustness"

********************************************************************************
* LOAD AND MERGE DATA
********************************************************************************

import delimited "fac_master_FINAL.csv", clear

gen yq = yq(year, quarter)
format yq %tq
tsset yq

gen fac_lag1 = L.score
label variable fac_lag1 "FAC score (t-1)"
label variable score "FAC score"

tempfile fac_data
save `fac_data'

use "fred_macro.dta", clear

rename fedfunds ffr
rename unemployment unemp

gen inflation_gap = inflation - 2
gen output_gap = -2 * (unemp - 5)

label variable ffr "Federal funds rate"
label variable inflation_gap "Inflation gap"
label variable output_gap "Output gap"

tempfile macro_data
save `macro_data'

use `fac_data', clear
merge 1:1 yq using `macro_data', keep(match) nogen

tsset yq
cap drop fac_lag1
gen fac_lag1 = L.score
label variable fac_lag1 "FAC score (t-1)"

********************************************************************************
* CREATE PHASE DUMMIES
********************************************************************************

* 24-month baseline (8 quarters)
local w = 8

gen precrisis = 0
replace precrisis = 1 if yq >= yq(1987, 4) - `w' & yq < yq(1987, 4)
replace precrisis = 1 if yq >= yq(1998, 3) - `w' & yq < yq(1998, 3)
replace precrisis = 1 if yq >= yq(2008, 3) - `w' & yq < yq(2008, 3)
replace precrisis = 1 if yq >= yq(2020, 1) - `w' & yq < yq(2020, 1)

gen crisis = 0
replace crisis = 1 if yq >= yq(1987, 4) & yq <= yq(1988, 1)
replace crisis = 1 if yq >= yq(1998, 3) & yq <= yq(1998, 4)
replace crisis = 1 if yq >= yq(2008, 3) & yq <= yq(2009, 2)
replace crisis = 1 if yq >= yq(2020, 1) & yq <= yq(2020, 2)

gen postcrisis = 0
replace postcrisis = 1 if yq > yq(1988, 1) & yq <= yq(1988, 1) + `w'
replace postcrisis = 1 if yq > yq(1998, 4) & yq <= yq(1998, 4) + `w'
replace postcrisis = 1 if yq > yq(2009, 2) & yq <= yq(2009, 2) + `w'
replace postcrisis = 1 if yq > yq(2020, 2) & yq <= yq(2020, 2) + `w'

gen normal = (precrisis == 0 & crisis == 0 & postcrisis == 0)

* 18-month windows (6 quarters)
local w = 6

gen precrisis_18 = 0
replace precrisis_18 = 1 if yq >= yq(1987, 4) - `w' & yq < yq(1987, 4)
replace precrisis_18 = 1 if yq >= yq(1998, 3) - `w' & yq < yq(1998, 3)
replace precrisis_18 = 1 if yq >= yq(2008, 3) - `w' & yq < yq(2008, 3)
replace precrisis_18 = 1 if yq >= yq(2020, 1) - `w' & yq < yq(2020, 1)

gen postcrisis_18 = 0
replace postcrisis_18 = 1 if yq > yq(1988, 1) & yq <= yq(1988, 1) + `w'
replace postcrisis_18 = 1 if yq > yq(1998, 4) & yq <= yq(1998, 4) + `w'
replace postcrisis_18 = 1 if yq > yq(2009, 2) & yq <= yq(2009, 2) + `w'
replace postcrisis_18 = 1 if yq > yq(2020, 2) & yq <= yq(2020, 2) + `w'

* 30-month windows (10 quarters)
local w = 10

gen precrisis_30 = 0
replace precrisis_30 = 1 if yq >= yq(1987, 4) - `w' & yq < yq(1987, 4)
replace precrisis_30 = 1 if yq >= yq(1998, 3) - `w' & yq < yq(1998, 3)
replace precrisis_30 = 1 if yq >= yq(2008, 3) - `w' & yq < yq(2008, 3)
replace precrisis_30 = 1 if yq >= yq(2020, 1) - `w' & yq < yq(2020, 1)

gen postcrisis_30 = 0
replace postcrisis_30 = 1 if yq > yq(1988, 1) & yq <= yq(1988, 1) + `w'
replace postcrisis_30 = 1 if yq > yq(1998, 4) & yq <= yq(1998, 4) + `w'
replace postcrisis_30 = 1 if yq > yq(2009, 2) & yq <= yq(2009, 2) + `w'
replace postcrisis_30 = 1 if yq > yq(2020, 2) & yq <= yq(2020, 2) + `w'

********************************************************************************
* TABLE B4: WINDOW LENGTH ROBUSTNESS (HC3 STANDARD ERRORS)
********************************************************************************

di ""
di "======================================================================"
di "TABLE B4: WINDOW LENGTH ROBUSTNESS (HC3 STANDARD ERRORS)"
di "======================================================================"

* Column 1: 18-month post-crisis
reg ffr inflation_gap output_gap fac_lag1 if postcrisis_18 == 1, vce(hc3)
outreg2 using "results/robustness/table_b4_window_robustness.doc", replace word ///
    title("Table B4: FAC Influence in Post-Crisis Periods by Window Length") ///
    ctitle("18-month") ///
    label dec(3) ///
    addtext(Phase, Post-crisis, Window, 18 months, SE type, HC3)

* Column 2: 24-month post-crisis (baseline)
reg ffr inflation_gap output_gap fac_lag1 if postcrisis == 1, vce(hc3)
outreg2 using "results/robustness/table_b4_window_robustness.doc", append word ///
    ctitle("24-month") ///
    label dec(3) ///
    addtext(Phase, Post-crisis, Window, 24 months, SE type, HC3)

* Column 3: 30-month post-crisis
reg ffr inflation_gap output_gap fac_lag1 if postcrisis_30 == 1, vce(hc3)
outreg2 using "results/robustness/table_b4_window_robustness.doc", append word ///
    ctitle("30-month") ///
    label dec(3) ///
    addtext(Phase, Post-crisis, Window, 30 months, SE type, HC3)

********************************************************************************
* TABLE B5: PHASE ANALYSIS (HC3 STANDARD ERRORS)
********************************************************************************

di ""
di "======================================================================"
di "TABLE B5: PHASE ANALYSIS (HC3 STANDARD ERRORS)"
di "======================================================================"

* Column 1: Normal periods (n=156, robust SEs acceptable)
reg ffr inflation_gap output_gap fac_lag1 if normal == 1, vce(hc3)
outreg2 using "results/robustness/table_b5_phase_analysis.doc", replace word ///
    title("Table B5: FAC Influence by Financial Cycle Phase") ///
    ctitle("Normal") ///
    label dec(3) ///
    addtext(Phase, Normal, Window, 24 months, SE type, HC3)

* Column 2: Pre-crisis
reg ffr inflation_gap output_gap fac_lag1 if precrisis == 1, vce(hc3)
outreg2 using "results/robustness/table_b5_phase_analysis.doc", append word ///
    ctitle("Pre-Crisis") ///
    label dec(3) ///
    addtext(Phase, Pre-crisis, Window, 24 months, SE type, HC3)

* Column 3: Crisis
reg ffr inflation_gap output_gap fac_lag1 if crisis == 1, vce(hc3)
outreg2 using "results/robustness/table_b5_phase_analysis.doc", append word ///
    ctitle("Crisis") ///
    label dec(3) ///
    addtext(Phase, Crisis, Window, 24 months, SE type, HC3)

* Column 4: Post-crisis
reg ffr inflation_gap output_gap fac_lag1 if postcrisis == 1, vce(hc3)
outreg2 using "results/robustness/table_b5_phase_analysis.doc", append word ///
    ctitle("Post-Crisis") ///
    label dec(3) ///
    addtext(Phase, Post-crisis, Window, 24 months, SE type, HC3)

********************************************************************************
* TABLE B6: BOOTSTRAP CONFIDENCE INTERVALS
********************************************************************************

di ""
di "======================================================================"
di "TABLE B6: BOOTSTRAP VALIDATION"
di "======================================================================"

* Store results for comparison table
tempname memhold
postfile `memhold' str20 phase n coef_ols se_hc3 ci_boot_lo ci_boot_hi using "bootstrap_results.dta", replace

* Normal periods
quietly reg ffr inflation_gap output_gap fac_lag1 if normal == 1, vce(hc3)
local n_normal = e(N)
local b_normal = _b[fac_lag1]
local se_normal = _se[fac_lag1]

quietly bootstrap _b[fac_lag1], reps(1000) seed(12345): reg ffr inflation_gap output_gap fac_lag1 if normal == 1
matrix ci = e(ci_percentile)
local ci_lo_normal = ci[1,1]
local ci_hi_normal = ci[2,1]

post `memhold' ("Normal") (`n_normal') (`b_normal') (`se_normal') (`ci_lo_normal') (`ci_hi_normal')

di "Normal (n=" `n_normal' "): coef = " %6.3f `b_normal' ", HC3 SE = " %6.3f `se_normal'
di "  Bootstrap 95% CI: [" %6.3f `ci_lo_normal' ", " %6.3f `ci_hi_normal' "]"

* Pre-crisis
quietly reg ffr inflation_gap output_gap fac_lag1 if precrisis == 1, vce(hc3)
local n_pre = e(N)
local b_pre = _b[fac_lag1]
local se_pre = _se[fac_lag1]

quietly bootstrap _b[fac_lag1], reps(1000) seed(12345): reg ffr inflation_gap output_gap fac_lag1 if precrisis == 1
matrix ci = e(ci_percentile)
local ci_lo_pre = ci[1,1]
local ci_hi_pre = ci[2,1]

post `memhold' ("Pre-crisis") (`n_pre') (`b_pre') (`se_pre') (`ci_lo_pre') (`ci_hi_pre')

di "Pre-crisis (n=" `n_pre' "): coef = " %6.3f `b_pre' ", HC3 SE = " %6.3f `se_pre'
di "  Bootstrap 95% CI: [" %6.3f `ci_lo_pre' ", " %6.3f `ci_hi_pre' "]"

* Crisis
quietly reg ffr inflation_gap output_gap fac_lag1 if crisis == 1, vce(hc3)
local n_crisis = e(N)
local b_crisis = _b[fac_lag1]
local se_crisis = _se[fac_lag1]

quietly bootstrap _b[fac_lag1], reps(1000) seed(12345): reg ffr inflation_gap output_gap fac_lag1 if crisis == 1
matrix ci = e(ci_percentile)
local ci_lo_crisis = ci[1,1]
local ci_hi_crisis = ci[2,1]

post `memhold' ("Crisis") (`n_crisis') (`b_crisis') (`se_crisis') (`ci_lo_crisis') (`ci_hi_crisis')

di "Crisis (n=" `n_crisis' "): coef = " %6.3f `b_crisis' ", HC3 SE = " %6.3f `se_crisis'
di "  Bootstrap 95% CI: [" %6.3f `ci_lo_crisis' ", " %6.3f `ci_hi_crisis' "]"

* Post-crisis (24-month)
quietly reg ffr inflation_gap output_gap fac_lag1 if postcrisis == 1, vce(hc3)
local n_post = e(N)
local b_post = _b[fac_lag1]
local se_post = _se[fac_lag1]

quietly bootstrap _b[fac_lag1], reps(1000) seed(12345): reg ffr inflation_gap output_gap fac_lag1 if postcrisis == 1
matrix ci = e(ci_percentile)
local ci_lo_post = ci[1,1]
local ci_hi_post = ci[2,1]

post `memhold' ("Post-crisis") (`n_post') (`b_post') (`se_post') (`ci_lo_post') (`ci_hi_post')

di "Post-crisis (n=" `n_post' "): coef = " %6.3f `b_post' ", HC3 SE = " %6.3f `se_post'
di "  Bootstrap 95% CI: [" %6.3f `ci_lo_post' ", " %6.3f `ci_hi_post' "]"

postclose `memhold'

* Save main dataset before loading bootstrap results
save "fac_window_robustness.dta", replace

* Create formatted Table B6
use "bootstrap_results.dta", clear
list, noobs clean

* Export to Word
gen ci_string = "[" + string(ci_boot_lo, "%6.3f") + ", " + string(ci_boot_hi, "%6.3f") + "]"
gen se_string = "(" + string(se_hc3, "%6.3f") + ")"

export delimited using "results/robustness/table_b6_bootstrap.csv", replace

********************************************************************************
* SUMMARY OUTPUT
********************************************************************************

* Reload main dataset (was replaced by bootstrap_results.dta above)
use "fac_window_robustness.dta", clear

di ""
di "======================================================================"
di "SUMMARY FOR PAPER"
di "======================================================================"
di ""
di "POST-CRISIS COEFFICIENTS BY WINDOW (HC3 SEs):"

quietly reg ffr inflation_gap output_gap fac_lag1 if postcrisis_18 == 1, vce(hc3)
local b18 = _b[fac_lag1]
local se18 = _se[fac_lag1]
local p18 = 2*ttail(e(df_r), abs(_b[fac_lag1]/_se[fac_lag1]))
local n18 = e(N)

quietly reg ffr inflation_gap output_gap fac_lag1 if postcrisis == 1, vce(hc3)
local b24 = _b[fac_lag1]
local se24 = _se[fac_lag1]
local p24 = 2*ttail(e(df_r), abs(_b[fac_lag1]/_se[fac_lag1]))
local n24 = e(N)

quietly reg ffr inflation_gap output_gap fac_lag1 if postcrisis_30 == 1, vce(hc3)
local b30 = _b[fac_lag1]
local se30 = _se[fac_lag1]
local p30 = 2*ttail(e(df_r), abs(_b[fac_lag1]/_se[fac_lag1]))
local n30 = e(N)

di "  18-month (n=" `n18' "): " %6.3f `b18' " (SE=" %5.3f `se18' ", p=" %5.3f `p18' ")"
di "  24-month (n=" `n24' "): " %6.3f `b24' " (SE=" %5.3f `se24' ", p=" %5.3f `p24' ")"
di "  30-month (n=" `n30' "): " %6.3f `b30' " (SE=" %5.3f `se30' ", p=" %5.3f `p30' ")"
di ""
di "PHASE COEFFICIENTS (24-month, HC3 SEs):"
di "  Normal (n=" `n_normal' "):     " %6.3f `b_normal' " (SE=" %5.3f `se_normal' ")"
di "  Pre-crisis (n=" `n_pre' "):  " %6.3f `b_pre' " (SE=" %5.3f `se_pre' ")"
di "  Crisis (n=" `n_crisis' "):      " %6.3f `b_crisis' " (SE=" %5.3f `se_crisis' ")"
di "  Post-crisis (n=" `n_post' "): " %6.3f `b_post' " (SE=" %5.3f `se_post' ")"
di ""
di "======================================================================"

save "fac_window_robustness.dta", replace

log close
